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ABSTRACT 

Magnetic reconnection in the partially ionized solar chromosphere is studied in 
2.5-dimensional magnetohydrodynamic simulations including radiative cooling 
and ambipolar diffusion. A Harris current sheet with and without a guide field 
is considered. Characteristic values of the parameters in the middle chromo¬ 
sphere imply a high magnetic Reynolds number of ~ 10®-10^ in the present 
simulations. Fast magnetic reconnection then develops as a consequence of the 
plasmoid instability without the need to invoke anomalous resistivity enhance¬ 
ments. Multiple levels of the instability are followed as it cascades to smaller 
scales, which approach the ion inertial length. The reconnection rate, normal¬ 
ized to the asymptotic values of magnetic field and Alfven velocity in the inflow 
region, reaches values in the range ~ 0.01-0.03 throughout the cascading plas¬ 
moid formation and for zero as well as for strong guide held. The out-flow velocity 
reaches ~ 40 kms“^. Slow-mode shocks extend from the X-points, heating the 
plasmoids up to ~ 8 x 10^ K. In the case of zero guide held, the inclusion of am¬ 
bipolar diffusion and radiative cooling both cause a rapid thinning of the current 
sheet (down to ~ 30 m) and early formation of secondary islands. Both of these 
processes have very little effect on the plasmoid instability for a strong guide 
held. The reconnection rates, temperature enhancements, and upward out-flow 
velocities from the vertical current sheet correspond well to their characteristic 
values in chromospheric jets. 
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Subject headings: magnetic reconnection - (magnetohydrodynamics) MHD - ra¬ 
diation: dynamics - Sun: activity - Sun: chromosphere 


1. Introduction 


It is well known that magnetic reconnection plays an important role in the dynamic 
solar corona, where it is a key process in eruptions (flares, filament eruptions, coronal mass 
ejections), jets, and other phenomena. In the lower solar atmosphere, magnetic reconnec¬ 
tion is also very common in different forms of activity, such as microffares (e.g., Brosius & 
Holman 2009; Gontikakis, Winebarger, & Patsourakos 2013), the flux cancellation process 
in the photosphere (e.g., Zwaanl987; van Ballegooijen & Martens 1989; Cameron, Vbgler, & 
Schussler2007), chromospheric jets (e.g., Chae, Moon, & Park 2003; Shibata et al. 2007; Liu 
et al. 2009; Shen et al. 2011; Shen et al. 2012 ; Morton 2012 ; Bharti, Hirzberger, & Solanki 
2013), and temperature enhancements in Ellerman bombs (e.g., Zachariadis, Alissandrakis, 
& Banosl987; Denkerl997; Dara et al.l997; Qiu et al.2000; Georgoulis et al. 2002 ; Fang 
et al. 2006). Magnetic reconnection probably plays an important role as the energy source 
for the heating of the chromosphere and corona (e.g., Parker 1972; Sturrock 1999; Klimchuk 
2006). The source of the mass and energy in the solar wind may also be related to magnetic 
reconnection in the lower solar atmosphere. 


Compared to the solar corona, the mass density and plasma fd in the chromosphere and 
photosphere are higher, and the temperature is much lower, mostly below 10^ K, except in the 
upper chromosphere. Therefore, the plasma in this height range is only partially ionized with 
an ionization degree of about 10“^-10“^ (e.g., Vernazza, Avrett, & Looser 1981; Pneuman, 


Solanki, & Stenflo 1986| Khomenko, Collados, & Felipe 2008). The magnetic diffusivity 
is of order 77 ~ ( 10 ^- 10 ^) m^s“^ in the range from the upper photosphere to the middle 
chromosphere, so that the Lundquist number (magnetic Reynolds number) for a length- 
scale L ~ 1 Mm and Alfven velocity va r\j ( 10 - 100 ) kms ^ is of order S = Lva/t} ~ 10 ®- 10 ®. 
The classical Sweet-Parker model for a current sheet of this length yields a reconnection 
rate ysp ~ ~ 10 “^- 10 “^, far smaller than the reconnection rates inferred for events 

of chromospheric activity. For example, Nishizuka et al. (2011) estimated 7 ~ 0.02-0.1 for 
chromospheric anemone jets from their life time and the estimated Alfven crossing time. 
The same estimate for type-II spicules, which have life times of 10-150 s, sizes of ~200 km, 
and ambient Alfven velocities based on the total density of ~ 100 kms“^ (de Pontieu et al. 
2007; Moore et al. 2011), yields 7 ~ 0 . 01 - 0 . 2 . Such and even larger disparity is typically 


found for weakly ionized astrophysical plasmas like, for example, the interstellar medium 
and protostellar and protoplanetary disks (Zweibel et al. 2011 and references therein). 
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Many numerical simulations applied to phenomena of solar activity employ ‘anomalous’ 
resistivity enhancements to obtain fast reconnection. However, the exact forms of anomalous 
resistivity used in magnetohydrodynamic (MHD) models are not directly deducible from 
the kinetic theory and simulations of real physical systems, and therefore some simplifying 
assumptions are still required for the model of anomalous resistivity (see, e.g., 
and references therein). 


Ni et al.||m2 


Steady-state reconnection rates in weakly ionized plasma can be significantly enhanced 
above the Sweet-Parker value by radiative cooling and by ambipolar diffusion; however, both 
are efficient only if the guide field is weak. Radiative cooling reduces the pressure in the 
current sheet, thus allowing the sheet to thin and, correspondingly, the current density to 
rise. This implies higher Lorentz forces accelerating the reconnection outflow, resulting in 
an enhancement over the Sweet-Parker rate of order where A is the compression factor 


of the sheet due to the cooling (Dorman & Kulsrud 1995; Uzdensky & McKinney 2011). In 


the incompressible strong guide held case, the cooling leads to a weaker enhancement of the 
reconnection rate through the higher resistivity, but in a chromospheric environment this can 
hardly reach an order of magnitude. Ambipolar diffusion in the current sheet decouples the 
ions and neutrals, so that only the ion pressure is available to balance the Lorentz force, which 
results in a strong current sheet thinning for weak ionization, allowing rapid flux annihilation 


(Brandenburg & Zweibel 1994). The enhancement over the Sweet-Parker reconnection rate 


can be very large in the special case of antiparallel fleld (Heitsch fc Zweib^|2003a), however. 


already a tiny guide fleld component suppresses this fast reconnection regime (Heitsch & 


Zweibel 2003b). 


In partially ionized plasma, fast reconnection also results from the recombination effect 
under non-equilibrium ionization conditions. The numerical results from a two fluid model 


in the papers by Leake et al. (2012, 2013) indicate that the reconnection rate can been 


increased to above 0.05 solely due to the decoupling of neutrals from ions. Therefore, when 
the neutral-ion collisional mean free path exceeds the current sheet width, the two-fluid 
model including the recombination effect should be applied in studying the reconnection 
process. 


Reconnection also becomes fast when the current sheet thins down to the ion inertial 
scale di where the Hall term is relevant (Mandt et al. 1994; Malyshkin & Zweibel |2011[ ). 
However, a dynamic regime of reconnection which involves plasmoids (magnetic islands) 
realizes high reconnection rates already at current sheet widths typically several orders of 
magnitude larger than di; this is now commonly referred to as the plasmoid instability 
(Loureiro et al. 2007; Daughton et al. 2009; Bhattacharjee et al. 2009). The motion of the 
plasmoids in the current sheet is largely constrained by ideal MHD (Finn & Kaw 1977); 
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hence, it is fast and largely independent of the resistivity, forming a highly elRcient internal 
engine of fast reconnection. The instability cascades to smaller scales such that current 
sheet sections between newly formed islands break up as well; the resulting high current 
densities in the increasingly thin current sheet fragments also facilitate a high reconnection 
rate. Moreover, this mechanism operates for current sheets with and without a guide field. 
Its high relevance for reconnection in the solar corona, especially for flares, has recently 
been demonstrated in a number of simulation studies (e.g., Murphy et al. 2012; Mei et al. 
2012 ; Guo et al. 2013). The plasmoid instability has recently been found to occur also in 


simulations of reconnection in weakly ionized plasma (Leake et al. 2012 ; Leake, Lukin, & 


Linton 2013[ ). Plasmoids indicating the action of the instability were detected in laboratory 
experiments (e.g., Dong et ^|2012 ), in the solar corona in the course of eruptions (e.g., Lin, 
Cranmer, & Farrugia 2008; Savage et al. 2010; Milligan et al. 2010; Takasao et al. 2012 ; 


Liu 2013), as well as in chromospheric jets (Shibata et al. 2007; Singh et al. 2012 ). 


The plasmoid instability occurs only if the Lundquist number and the aspect ratio L/5 
of the current sheet are sufficiently large, where 5 is the current sheet width. Both are related 
through the Sweet-Parker scaling for the current sheet width, (jsp ~ For fully ionized 

plasma, the critical Lundquist number, Scr, typically lies in the range of several 10^ to several 
10^, decreasing for increasing plasma /5 in the inflow region (e.g., Bhattacharjee et al. 2009; 
Huang & Bhattacharjee 2010; Ni et al. 2010; Ni et al. 2012; Ni et al. 2013). The minimum 
aspect ratio correspondingly is of order ~ 10^. Fast reconnection at a rate 7 ~ Scr^^'^ was 
found for S > Scr- Leake et al. ( 2012 , 2013) And the onset of the plasmoid instability in 


partially ionized plasma in basic agreement with this value. 


Since the characteristic Lundquist number estimated for chromospheric reconnection 
events is much higher than the critical Lundquist number for onset of the plasmoid instability, 
and since this process yields high reconnection rates for weak and strong guide fields, we 
focus on the plasmoid instability in chromospheric reconnection in the present investigation, 
considering both the regimes of strong and of vanishing guide held. Parameters characteristic 
of the middle solar chromosphere and a low plasma P in the inflow region are considered. 
The energetically and dynamically important processes of radiative cooling and ambipolar 
diffusion are included, and their effect on the reconnection rate is compared with the effect 
of the cascading plasmoid instability. The formation of slow-mode shocks associated with 
the X-points between the plasmoids is demonstrated. Finally, we compare the obtained 
reconnection rates, upward reconnection outflow velocities, and temperature enhancements 
with the values typically seen in or inferred for events of chromospheric activity. 


The model and numerical method are described in the following section. In Section 3, we 
present our numerical results and compare them with observations of chromospheric activity 
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phenomena. A summary and discussion are given in Section 4. 


2. Model and Numerical Method 


We only consider the hydrogen gas, so that the plasma is composed of neutral hydrogen 
atoms, ions (protons) and electrons. In general this system is described by the three-fluid 
MHD equations. However, a single-fluid approximation adding up all three equations can be 
used if the collisional coupling between ions and neutrals is strong, which is valid on length- 
scales exceeding the neutral-ion collisional mean free path Ani. The relevant length-scale here 
is the current sheet width S. We will find below that parameter values characteristic of the 
middle chromosphere yield critical current sheet widths for onset of the plasmoid instability 
several orders of magnitude above Ani, indicating that the onset and initial evolution of the 
instability in this part of the solar atmosphere can generally be described using the one- 
fluid equations. Since the instability tends to cascade to small scales comparable to the ion 
inertial length, which can be of similar magnitude as Ani or smaller (depending on the ion 
density), it is clear that a complete description of plasmoid-mediated reconnection in the 
chromosphere will require a multi-fluid treatment (such as recently developed by 

. However, the present investigation demonstrates that the plasmoid instability 
yields high reconnection rates already within the validity range of the one-fluid 
approximation adopted here. Since the neutrals are coupled to the ions in this approximation, 
the kinematic behaviour of the plasma bears a considerable degree of similarity to the fully 
ionized case, but radiative losses and ambipolar diffusion can nevertheless influence the 
dynamics strongly. 


2012, 2013 


occurs and 


Leake et al. 


Additionally, we approximate the gas pressure tensor as a scalar, assume isothermal 
conditions (for simplicity and in the absence of detailed information to the opposite), 
Ti = Tn = Tg, and neglect heat conduction. The latter is based on the fact that heat 
conduction is energetically less important than radiative cooling for the high densities and 
low temperatures characteristic of the middle chromosphere (Lin et al. 1992). Additionally, 
heat conduction into and out of plasmoids is largely across the magnetic field, especially 
in a two-dimensional description, thus, it tends to be small (see Section for a detailed 
discussion). This results in the set of basic equations (see, e.g., Khomenko &: Collados||2012) 
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As usual, p is the plasma mass density, v is the centre of mass velocity, e is the total energy 
density, B is the magnetic field, p is the magnetic diffusivity, and p is plasma thermal 
pressure, g = —273.9 ms“^ By is the gravitational acceleration of the Sun. The ambipolar 
electric field is given by Eau = Pq^Pad[C^ x B) x B] x B, where pad is the ambipolar 
diffusion coefficient. £rad is the radiative cooling function and PL is the heating function, 
li = Ue/uH is the ionization degree of the plasma. In the above equations, the following 
definitions are used: 
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( 10 ) 


The subscripts ‘n, i, e, H’ refer to neutral hydrogen atoms, ions, electrons and the total num- 
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ber of hydrogen particles (neutral and ionized), respectively, such that the number densities 
are related by uh = Un + Ui = Un + Ug. 


In computing the radiative loss we follow Gan & Fang (1990) and set 


Aad = -1.547 X 10-"^ 1^0 ( ^ 


P 


( 11 ) 


where the subscript ‘0’ represents the initial value at time t = 0. The coefficient a in 


the above function is the same as the one in Equation (11) in Gan & Fang (1990). Their 
model is applicable up to T ~ 10® K and is well supported by recent investigations of the 
active solar chromosphere (e.g., Fang, Chen, & Ding 2003; Jiang, Fang, & Chen 2010; Xu 
et al. 2011). We assume the plasma to be in ionization equilibrium. The cooling processes 
considered in this model ( Raymond et al.||1976 ) are permitted, forbidden and semiforbidden 
line transitions, including contributions from dielectric recombination and bremsstrahlung, 
radiative recombination, and two-photon continua. The heating function is chosen to initially 
balance the radiative loss exactly. 


H = 


1.547 X 10-^‘^Yio 


m\ 


aTr 


1.5 


( 12 ) 


The simulation domain extends from a: = 0 to a: = Lq in ^ direction and from y = 0 to 
y = 2Lo in y direction, with Lq = 10® m. Open boundary conditions are used in both x and 
y directions. 

Three models of the initial current sheet are simulated in this work, differing in the guide 
field strength and in the inclusion of a density stratification and gravity. For each model 
we consider three cases. Cases A-C, D-F, and G-I in Models I, II, and II, respectively. In 
Cases A, D and G radiative cooling, heating, and ambipolar diffusion are excluded. Cases B, 
E and H include heating and radiative cooling but not ambipolar diffusion. Cases C, F and 
I include ambipolar diffusion but not the heating and radiative cooling terms. The other 
initial and boundary conditions and the physical parameters are identical within the cases 
for each model. 


The computations are performed using the MHD simulation code NIRVANA (version 
3.6; Ziegler 2011). This code does not yet have the capability to describe ionization and 
recombination, so we work here with a non-time dependent ionization ratio Vo and defer the 
inclusion of these effects to a follow-up investigation. 


Adaptive mesh refinement (AMR) is applied. We start the simulation from a base- 
level grid of 160 x 320. The highest refinement level is 10, which corresponds to a grid 
resolution /S.x ~ 6.1 m. This resolution is comparable to the neutral-ion collisional mean 












fee path Ani in the middle solar chromosphere (see Table . Convergence studies have been 
carried out by repeating some simulations with both a lower and a higher resolution, with 
the highest rehnement level respectively limited to 9 and 11. The numerical results in the 
lower resolution case and the higher one are very similar. 


2.1. Models I and II 


Since the effects of radiative cooling and ambipolar diffusion on the reconnection rate 
depend strongly on the guide field ( Uzdensky h McKinney||2011 ; Heitsch fc Zweibe'I||2003b ), 
we consider the Harris sheet with strong guide field (Model I) and without a guide field 
(Model II). These are realized by employing two versions of the Harris current sheet equi¬ 
librium, both negelecting gravity and oriented vertically. A force-free Harris current sheet, 
which is appropriate in a low-/5 environment and has a strong guide field by its nature, is 
used as Model I, 


S.0 = 0 (13) 

Byo = 6otanh[(a: — Lo/2)/0.05Lo] (14) 

B,o = 6o/cosh[(x-Lo/2)/0.05Lo], (15) 

where bo = 0.01 T. The initial current sheet width thus is (5o = O.ILq. Due to the force- 

freeness and neglect of gravity, the initial equilibrium thermal pressure is uniform. The 
initial plasma (5 is also uniform and set to /3o = 0.1. This yields the initial plasma thermal 
pressure po = 12.57r“^ Pa. 


To realize a model without guide field, we must deviate from force freeness. Model H 
uses a standard Harris current sheet in equilibrium with a pressure gradient, 

5x0 = 0 (16) 

Byo = 6otanh[(x — Lo/2)/0.05Lo] (17) 

5,0 = 0. (18) 


From equation (2), the initial plasma pressure in Model H is 

Po = (1+ ^o - (tanh[(x - Lo/2)/0.05Lo])^);^, 

zno 

and we set the initial asymptotic plasma (5 also to /5o(|x| oo) = 0.1. 
In all cases we assume the initial equilibrium temperature as 

To = Toft [1 -|- (1 -|- tanh(12(p — To)/To))/2], 


( 19 ) 


( 20 ) 
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where To;, = 3500 K. The functions Pq{x) and To(r/) yield the total mass density pQ{x,y) 
and particle number density uho + for each case. Then we calculate the initial ionization 
degree = '^eo/'^no using equation (12) in Gan & Fang (1990). Our settings yield rather 
similar values of the initial gas pressure, mass density, held strength, and Alfven velocity 
uao = in the inhow region (x = 0 and Lq, and y = Lq) for both models. 

However, in Model II the gas pressure and mass density in the current sheet region are 
higher by a factor 11 compared to the inhow region. The resulting Alfven velocity in the 
current sheet is lower than in Model I by a factor 3. We should point out that the Alfven 
velocity calculated in this work is based on the total plasma density, as is appropriate for 
strongly coupled partially ionized plasmas (at scales S ^ Ani); see, e.g., Zweibel (1989). 


The basic plasma parameters in the current sheet region are compiled for the two models 
in Table (along with the values of the resulting Lundquist number based on Lq, Alfven 
velocity, and length-scales); they correspond to the conditions in the middle chromosphere, 
i.e., at roughly 1 Mm above the photosphere, with the field strength being representative 
of active regions. The neutral-ion collisional mean free path is given by the ratio of the 
thermal velocity of the neutrals and the neutral-ion collision frequency, Ani = with 

VTn = z/ki = niCTki a/ 8A:BTki/(7rmki), = (Tk -F Ti)/2, mu = mkmi/(mk + mi), 

(where the subscripts k and 1 denote the species), and the value of the collisional cross section 
(Jki is taken from Khomenko fc Collados| (|2012). 


Initial perturbations of both magnetic field (B', with V • B' = 0) and velocity (v') are 
applied at t = 0 in all cases to trigger the reconnection process. The perturbations result in 
a thinning of the current sheet in two sections between a set of three primary islands, whose 
midpoints are located at y = 0, Lq, and 2Lo (see Fig. 1(a)). In this paper, we will focus 
only on the section in the domain Lq < y < 2Lq, i.e., the bottom half of the box is used 
as an auxiliary part of the computation only. Its function is to yield a stationary primary 
plasmoid at the bottom of the height range of interest, which is not influenced by any effects 
of a numerical boundary; the primary plasmoid acts like a line-tied bottom of the current 
sheet of interest. The initial temperature Tq varies from 5250 K to 7000 K in the domain of 
interest, but it quickly increases to the upper value above the stationary plasmoid at y = Lq, 
such that essentially all of the upper current sheet initially is at Tq ^ 7000 K. The resulting 
ionization degree in the current sheet is 0.5% in Model I and 0.2% in Model II. 


We use a simple parametrization of the magnetic diffusion which relatively closely 
matches the diffusion computed from a model atmosphere in Khomenko & Collados (2012), 
r) = b X 10^(3500/T)^'^ m^s“^. This yields the initial diffusion in the current sheet as 
1.77 X 10^ m^ s The Lundquist number based on this diffusivity, uao, and Lq, which corre¬ 
sponds to the ‘global scale’ of the current sheet in the upper part of the box, is Sq = 1.88 x 10® 
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in Model I and 5.66 x 10® in Model II. rj decreases as the temperature increases. From the 
numerical results presented in the following section, the value of rj at the main X point drops 
to ~ 2000 m^ s“^ during the secondary instability process, and the Lundquist number then 
reaches lOF 


The ambipolar diffusion coefficient t^ad is defined as (Braginskii 1965) 

^AD ^ ^ ^ (Pi^in T Pe^en) ) 


( 21 ) 


where, respectively, r'in and z^en are the ion-neutral and electron-neutral collision frequencies 
defined analogous to above, with the cross sections (Tin = 5 x 10“^® m^ and cxen = 10“^® m^ 


again taken from Khomenko & Collados (2012). This yields 


Pad = 1.65 X 10 - 1)-^^ 


m® skg 


,-i 


( 22 ) 


with the temperature inserted in Kelvin. The ambipolar diffusion coefficient is a function of 
ionization degree, temperature and plasma density. According to the model atmospheres by 


Khomenko, Collados, & Felipe (2008) and Pneuman, Solanki, & Stenflo (1986), the magni¬ 
tude of Pad varies from 10 to 10® m® skg“^ from the bottom to the top of the chromosphere. 


Here we set pad = 5 x 10^(ro/T)^/^(po/p)^ m®skg 
diffusion coefficient in the middle chromosphere. 


This is representative of the ambipolar 


For the radiative loss and heating function in Models I and II, we simplify the coefficient 
a compared to the expression in Gan fc Fang| ( |1990| ) , which is slowly varying in the middle 
chromosphere (about 1000 km above the photosphere), the height range of interest in the 
present study. Since we do not investigate the height dependence of the reconnection rate, we 
here use a representative uniform value of a = 0.01, which is chosen such that the resulting 
magnitude of the initial radiative loss is of order 0.07-3.5 Jm“®s 
results. 


® ^ consistent with their 


2.2. Model III 

Since the active section of the current sheet in our computations, Lq < y < 2Lo, is of 
order 1 Mm, it encompasses a considerable range of density and, correspondingly, Alfven 
velocity in the chromosphere. In Model III we therefore include gravity and the resulting 
stratification of thermal pressure and plasma density. It is not possible to construct an 
equilibrium consisting of a vertical Harris sheet without a guide field in the presence of 
gravity. Therefore, Model HI combines the force-free Harris sheet of Model I with density 
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stratification, gravity, and a balancing pressure gradient, given by 


dxPo{y) = 0, 

(23) 

dyPoiy) = -273.9po(l/), 

(24) 

(1 + ^0)PH0(l/) rp 

Po{y) — 

rrii 

(25) 


where, for simplicity, the initial temperature Tq = 7000 K is chosen to be uniform, and 
Pho = nnorrii. According to Equation (12) in Gan & Fang (1990), the ionization degree is 
simplified as 


Do = 


^eO 

^HO 


00 


(26) 


where 0o — 1-21 x 10^° m ^ and 0o ^ho in our model. From Equations (22-25), we obtain 


nm ^ ^Hoo exp 


4.706y\ 

■ Ao j’ 


(27) 


where we assume uhoo = m“'^. The same perturbations as those in Models I and II are 

applied; we also only focus on the section in the domain Lq < y < 2 Lq. From Equations (24- 
26), we calculate the initial distributions of hydrogen number density uho, ion number density 
Uio, ionization degree Alfven velocity uao, ion inertial length dio, and neutral-ion collision 
mean free path Anio- The values of these parameters in Model III at y = Lq, y = 1.5Z/o and 
y = 2Lo are compiled in Table [T] 

The initial Lundquist number So = Luao/p is of order 10® (see Table [^. The initial 
ambipolar diffusion coefficient takes values of tjad — 75.67 m^skg“^ at y = Lq, 2.57 x 
10^ m^skg“^ at y = 1.5Lo, and 8.71 x 10^ m^skg“^ at y = 2.0Lo- 


Since the density stratification in y direction is included in Model III, we here use the 
full y-dependent expression for a from Gan & Fang (1990), so that the resulting height 
dependence of the radiative loss is similar to their results. Gompared to their expression, 
we here increase the value of a by a factor of 10, such that the initial radiative loss is in 
the range ~ (0.2-1.3) Jm“^s“^ for Lq < y < 2Lo, considerably higher than the value of 
in Model I. This choice is made to demonstrate that a main result we will 


-3 ^-1 


~ 0.07 J m 

obtain for Model I, namely that radiative losses in the chromosphere do not significantly 
influence the evolution of current sheets with a strong guide field, remains valid for the much 
higher level of radiative losses. The levels chosen for the initial radiative losses in Models I 
and III both fall in the range deduced in Gan & Fang (1990). 


From the above descriptions, we find that model III is close to the conditions in the 
middle chromosphere, i.e., in the height range ~ (0.5-1.4) Mm above the photosphere. 
Figure 1(b) displays the initial magnetic field and plasma density. 
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3. Results 


3.1. Cascading plasmoid instability with internal slow-mode shocks 


The Lundquist number in our models exceeds the critical Lundquist number for onset 
of the plasmoid instability substantially. As the Harris sheet evolves to a thinner current 
sheet in response to the initial perturbation, it first passes through a Sweet-Parker-like phase 
(with the exception of Case F) as indicated by its smooth, very elongated structure (which 
contains an X-point at y ~ 1.5 Lq) and an average reconnection rate 7 ~ 10“^ in agreement 
with the Sweet-Parker scaling (see below for the method of quantitative analysis). When 
the current sheet sufficiently thins such that a critical aspect ratio is reached, the plasmoid 
instability sets in, forming multiple islands. Subsequently, the instability cascades to smaller 
scales. Throughout this process a high reconnection rate, rapidly fluctuating in the range 


~ 0.01-0.03, is realized (see Sect. 3 . 2 ). The plasmoid cascading process is found to be very 
similar in all nine cases run. 

Fig. 2 displays the plasmoid cascading process for Model III, Case I with ambipolar 
diffusion included. The left panel shows the primary island at y = Lq produced by the 
initial perturbation and a number of secondary islands formed by the plasmoid instability 
(note that the primary island at y — 2 Lq was ejected by the slow reconnection outflows 
that evolved in the Sweet-Parker-like phase). As we zoom in (panels 2-3), several third- 
order plasmoids become visible in the current sheet, which has further narrowed between the 
secondary plasmoids. Zooming in further (panels 3-4), fourth-order plasmoids can be seen 


clearly. This plasmoid cascading process is very similar to the results of Barta et al. (2011), 


who have simulated it for coronal parameters. It is found here for the first time in partially 
ionized plasma. 

The half-width of the thinnest fragment current sheet in Fig. 2 is about 25 m, clearly 
larger than the neutral-ion collisional mean free path Ani (Table [^, thns consistent with 
the one-fluid approximation. The cascading plasmoid instability occurs in all nine runs 
considered in this paper, reaching comparable reconnection rates in the cases with and 
without a guide field, which differs strongly from the results for stationary reconnection 


in partially ionized plasma (Heitsch & Zweibel 2003b; Uzdensky & McKinney 2011). The 


gravity of the Sun also has only little effect on the time dependent reconnection rate (see 
below). 

Given the large body of experience from one-fluid, Hall-MHD and kinetic models of 


reconnection in the plasmoid-unstable range (e.g., Huang, Bhattacharjee, & Sullivan 2011 
iDaughton et al. 2009), one can reasonably expect that the plasmoid cascading process 


tends to continue to even smaller scales than resolved here, down to the ion inertial length di. 
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In partially ionized plasma the process is modified by recombination within the plasmoids, 
reducing their current and fiux ( Leake et al.j2012 ), and it is not yet known how this infiuences 
the cascading process. A less dynamical evolution of the plasmoids is conceivable, but 
since the cascading occurs between previously formed plasmoids, it appears unlikely that 
recombination within plasmoids would suppress it completely. This is indeed indicated by 


the results in Leake et al. (2013), who find that that the plasmoid instability can raise the 


reconnection rate also beyond the validity range of the one-fluid approximation. We thus 
expect that the multiple levels of the plasmoid instability process connect the global, large 
MHD-scale reconnection with the local, small kinetic-scale reconnection also in partially 
ionized plasma. It is well known that the plasmoid instability in fully ionized plasma yields 
a high reconnection rate throughout the cascading process, and our simulations demonstrate 
that the cascading plasmoid instability operates in weakly ionized plasma in the one-fluid 
limit in basically the same manner as in fully ionized plasma, reaching a high reconnection 
rate as well. 

As is well known, reconnection at A-points always involves slow-mode shocks in the 
Petschek regime. In our simulations, we find that many small-scale slow-mode shocks form 
transiently at the edge of secondary plasmoids, attached to the neighboring A-point in a 
secondary fragment of the current sheet. Fig. 3 shows a pair of such shocks. Sudden changes 
of the magnetic field direction and current density at the left and right edges of the plasmoid 
can clearly be seen; these are very similar to the slow-mode shock structure displayed, e.g., 
in Forbes (1988), Schumacher & Kliem (1997) and Mei et al. (2012). Fig. 3(b) shows that 


the temperature increases strongly not only within the current sheet fragment but also 
downstream of the slow-mode shocks in the plasmoid. The angle between the shock front 
and the y-direction is ~ 5.4°. Fig. 4 displays the magnetic field and the current density 
along a cut in the x-direction at y = 1.691 Lq. One can see that the magnitude of the field 
component tangential to the shock, B\\, decreases rapidly toward the downstream side and 
that the current density A has a peak where changes rapidly, but the component normal 
to the shock, B±^ stays nearly uniform along the cut. This is exactly the behavior of a 
slow-mode shock. Similar structures exist in the outflow region of the current sheet. These 
slow-mode shocks could be a candidate mechanism to explain the chromospheric heating 


(Kumar, Kumar, & Uddin 2011). 
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3.2. Effects of radiative cooling and ambipolar diffusion 

3.3.L Strong guide field 

We now compare the simulation results of the three cases which were run for Model 1. 
The reconnection rate is computed as the rate of change of the magnetic flux accumulated 
between the 0 -point in the primary island at y = Lq and the main reconnection X-point 
(see Ni et al. 2012; Ni et al. 2013), 

. X ^ diipxit) - i^ojt)) 1 

dt boVAo' 

Here the flux function xp is defined through the relations = —dip/dy, By = dxp/dx, and 
the main reconnection X-point is determined as the X-point which has the highest tp value 
of all X-points in the box. This analysis is performed on the relatively low grid refinement 
level 3 which ensures uniform coverage of the current sheet section from the primary 0-point 
to the upper boundary throughout the runs. This choice still captures much of the temporal 
variability of the reconnection rate. Local quantities like the current density, temperature, 
and current sheet width at the main reconnection X-point are determined at the highest 
refinement level chosen by the code, to evaluate their full dynamic range; after the onset 
of the plasmoid instability this is level 10 in all nine cases. The current sheet width is 
determined as the full width at half maximum (FWHM) of the current density profile in x 
direction. 

The reconnection rate for Cases A-C is displayed in Fig. 5(a). The three cases are 
very similar to each other and show a dynamical behaviour that is basically analogous to 
reconnection in fully ionized plasma, as expected from our one-fluid approximation, especially 
for Case A in the absence of heating, radiative cooling, and ambipolar diffusion. Up to 
t ~ 26 s the reconnection rate remains at a low level, reaching the Sweet-Parker value of 
7 ~ ^ 10“^ for f > 20 s. At this time the decreasing current sheet width has also 

reached the Sweet-Parker value 5x ~ ~ 10^ m. This appears to be a Sweet-Parker 

reconnection phase, which is also supported by the elongated, smooth appearance of the 
current sheet (see Fig. [^. 

Subsequently, a fast rise of the reconnection rate commences, which quickly saturates 
at a level 7 ~ 0.02. The first secondary islands due to the plasmoid instability have clearly 
developed by the end of the fast rise, which we thus interpret as the linear phase of the 
plasmoid instability. All quantities show the high variability characteristic of the plasmoid 
instability in the further evolution. The current sheet width, measured at the main X-point, 
gradually decreases in the Sweet-Parker phase and enters a rapid decrease down to ha: ~ 30 m 
simultaneously with the onset of the plasmoid instability (Fig. 6 (a)). This minimum current 
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sheet width is comparable to the grid scale at rehnement level 10 {5x ~ 5Ax) but still much 
larger that the neutral-ion collisional mean free path and the ion inertial length (Table 1), so 
that our one-fluid approximation is valid throughout the computation. The current density 
at the main X-point shows a rather close inverse proportionality to the current sheet width. 
The average reconnection rate remains at a high level throughout the cascading phase of 
the instability. It shows a trend of gradual decrease to an average level of ~ 0.015 which is 
most likely due to the shortening of the active section of the current sheet as the primary 
island at y = Lq grows; the held lines in the inhow region must then bend increasingly before 
reconnecting at the X-point. The temperature at the main X-point increases rapidly from 
the onset of the plasmoid instability, saturating at ~ 3 x 10^ K (Fig. 7(a)). The corresponding 
decrease of the magnetic diffusivity raises the effective Lundquist number in the course of 
the plasmoid instability by nearly an order of magnitude. It should be noted that the peak 
temperature in the current sheet, located within the islands downstream of the slow-mode 
shocks, continues to rise, reaching ~ 8 x 10^ K during the plasmoid instability. 


The comparison of Cases A and B shows that the radiative cooling hardly has any effect 
on the evolution of the reconnecting current sheet at the chosen parameters of the middle 
chromosphere if the sheet includes a considerable guide field. This is fully in agreement 
with earlier investigations (e.g., Uzdensky fc McKinney||2011 ), which had pointed out that a 
sufficiently strong guide field suppresses the current sheet thinning due to radiative cooling 
by making the sheet incompressible. The only remaining effect on the reconnection rate is 
given by the temperature dependence of the magnetic diffusivity. [Uzdensky fc McKinney 


(2011) considered the case that radiative cooling in a stationary Sweet-Parker sheet with 
negligible heating increases the reconnection rate in this way. However, heating is a key 
factor of the energy balance in the chromosphere, and so our model also includes a heating 
term which is adjusted to match the initial cooling rate and to have the same dependence 
on the density. Consequently, the radiative losses will not cool the plasma below the initial 
temperature Tq; they have a strong effect only when the temperature increases considerably 
above Tq. Lower temperatures can be reached by adiabatic cooling if the plasma expands 
(an effect seen below in Fig. |^, but this is not the case here at the X-point, where the 
temperature has the same value as in Case A during the Sweet-Parker phase (Fig. 7(a)). The 
subsequent fast reconnection dominated by the plasmoid instability is largely independent 
of the magnetic diffusivity (Huang & Bhattacharjee 2010). Temperatures of (4-5)ro are 
reached in this phase, as in Case A, which shows that the strong heating at the slow-mode 
shocks dominates the radiative losses in this phase. 


Considering the effect of ambipolar diffusion in Case C, again only a weak influence on 
the reconnection rate is found. The only differences to Cases A and B are a more impulsive 
and slightly earlier onset of the plasmoid instability and a somewhat higher variability near 
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the end of the Sweet-Parker phase. We interpret this behaviour as follows. In the Sweet- 


Parker phase a strong thinning of the current sheet is inhibited by the guide field (Heitsch & 


Zweibelj2003b ). The ambipolar electric field nevertheless increases strongly as the sheet thins 
in response to the initial perturbation. However, the area of enhanced £1 ad does not reach the 
plane x = 0.5 Lq where the antiparallel field component reverses sign, but rather forms two 
layers on either side of of the field reversal plane (Fig. 0: a), (c)). Thus, it does not amplify 
the annihilation of the field in this plane. The width of the layers, ~ 500 m, agrees with the 
ion-neutral decoupling length scale Lad = where is the Alfven velocity based 

on the ion density (Zweibel et al. 2011). As a dominant A-point develops, the associated 
reconnection flows convect the two layers with enhanced ambipolar electric field to the field 
reversal plane (Fig. Bb),( c)), so that it now contributes to the change of flux (Equation]^. 
The comparison of Figs. [II|(a) and (b) shows that the rate of change of the magnetic field in 
the Sweet-Parker-like current sheet is dominated by the resistive contribution to the electric 
field, whereas it is dominated by the ambipolar electric field during the rapid rise of the 
reconnection rate at the first A-point of the commencing plasmoid instability. The layers 
of enhanced ambipolar field intermittently approach the field reversal plane already during 
the considerable thinning of the current sheet in the final part of the Sweet-Parker phase 
(Fig. 6(a)), causing the intermittent moderate enhancements of the reconnection rate seen 
in Fig. 5(a). Since the inclusion of the ambipolar diffusion terms causes a strong decrease 
of the adaptive time step, this run was terminated after the reconnection rate reached the 
saturation level. 


The three cases also show a very similar current sheet aspect ratio at the onset of the 
plasmoid instability. We refer to this quantity as the “onset aspect ratio” and determine it 
from the current sheet width 5x and the similarly defined FWHM length of the current sheet 
in y direction, 5y. We find 5x ~ 400 m and 5y ~ 3.5 x 10® m, i.e. an onset aspect ratio 
fions = SyjSx ~ 875 in all three cases (see Fig. 8 for Case B). This value lies considerably above 
the critical aspect ratio for onset of the plasmoid instability at the point of transition between 
the Sweet-Parker and plasmoid-mediated reconnection regimes, Scr , when the Lundquist 

10®-10^ (Huang & Bhattacharjee 2010; Ni 


number is varied to reach the critical value Scr ^ 


et al. 2012; Leake et al. 2012). Our simulations are performed with far higher Lundquist 
numbers in the range S ~ 10®-10^ (varying in time due to the evolving Sx and T), and 
correspondingly the current sheet aspect ratio rises considerably above the critical value 
Slr'^, roughly aproaching 


In summary, in current sheets with a strong guide field, the plasmoid instability mediates 
a fast reconnection regime in partially ionized plasma, while radiative cooling and ambipolar 
diffusion have only very little effect. For the parameters of the middle chromosphere, the 
reconnection rate is enhanced by a factor ~ 20 above the Sweet-Parker rate, reaching the 



















range of observationally inferred values. 


3.2.2. Zero guide field 

Next we consider the Harris current sheet with vanishing guide field (Model II, Cases D- 
F). A major difference to all corresponding runs for Model I is the much slower evolution 
of the current sheet, due to our choice to realize the necessary pressure enhancement in 
the sheet by a density enhancement and the resulting smaller Alfven velocity in the sheet. 
The reconnection rate, which is normalized to the field strength and Alfven velocity in the 
inflow region, thus takes a much longer time to rise. The density in the current sheet 
gradually decreases during our runs, and eventually the Alfven velocity in the sheet becomes 
comparable to the external Alfven velocity. The reconnection rate then reaches similar values 
as in Model I. 

In all three cases run for Model II, reconnection rates exceeding the Sweet-Parker value 
develop only from the onset of the plasmoid instability (Fig. |^, which happens at t ~ 92, 
84, and 50 s for Cases D, E, and F, respectively (Figs. 5(b) and 8 ). A Sweet-Parker-like 
regime develops only transiently, for ~ 10 s prior to the onset of the plasmoid instability, in 
Cases D and E as the current sheet width approaches the Sweet-Parker value (Fig. 6 (b)). As 
for Model I, the onset of the plasmoid instability is associated with a rapid decrease of the 
current sheet width. Due to the high inertia of the plasmoids in the initially dense current 
sheet, the instability shows a far more gradual overall development compared to the strong 
guide field cases. 

The reconnection rate in Case D without radiative cooling and ambipolar diffusion 
does not reach saturation but rather rises essentially monotonically (apart form the small 
time-scale fluctuations characteristic of the instability) until the main X-point leaves the 
box through the upper boundary at t ~ 230 s. A peak reconnection rate of 7 ~ 0.025 is 
then reached; it would be somewhat higher (lower) if the box were chosen somewhat larger 
(smaller). The current sheet width continues to decrease in the course of the cascading 
plasmoid instability until it reaches 5x ~ 30 m, as in Model I. The temperature at the 
A-point reached at the end of the run, T ~ 3 x 10^ K, is also similar to the result for 
Model I. 

The evolution in Case E with radiative cooling and heating is qualitatively similar to 
Case D, but it commences earlier and is more impulsive. Here the current sheet width drops 
to 5x ~ 30 m, and the reconnection rate rises to 7 ~ 0.025, in about one half of the time. 
This more dynamic behaviour is related to a faster decrease of the density in the current 
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sheet after t ~ 100 s. It appears that the reconnection rate saturates at this level until it 
sharply drops as the main X-point leaves the box at t ~ 170 s. The effect of the radiative 
cooling again remains minor prior to the onset of the plasmoid instability. The current sheet 
is then only slightly thinner than in Case D (Figs. 6 (b) and 8 ). The effect of radiative cooling 
becomes more obvious at the high temperatures in the course of the plasmoid instability, 
limiting the temperature to about half the peak value found in Case D (Fig. 7 (b)). Radiative 
cooling in Case D acts more efficiently than in Case B because here the densities are higher. 


As expected in the absence of a guide field, the inclusion of ambipolar diffusion in 
Case F allows the current sheet to thin very efficiently (Brandenburg & Zweibel 1994). The 
evolution here begins with a rapid decrease of the current sheet width at the main X-point 
down to (ix ~ 30 m (within ~ 10 s), which is near its smallest value observed in all runs at 
the AMR refinement levels 10 and 11. A thinning by a factor 500 is expected from the ratio 
of ion and neutral pressure (the inverse ionization degree for our isothermal conditions) when 
only the ion pressure in the sheet balances the magnetic pressure of the external region. The 
actual thinning goes beyond that by a factor ~ 3, but this already includes the onset of the 
plasmoid instability at t ~ 50 s and 5x ~ 80 m. The instability seamlessly grows out of the 
current sheet collapse mediated by the ambipolar diffusion, without leaving room for a Sweet- 
Parker-like phase, since the onset aspect ratio is obviously reached during the collapse. The 
reconnection rate grows in a manner similar to Cases D and E, at least to a level 7 ~ 0 . 02 . 
It is possible that it would grow to higher values if the run could be continued. However, 
the main reconnection X-point has moved out of our simulation domain after t = 150 s. It 
is seen that the reconnection rate then suddenly decreases to a low value, similar to Cases E 
and D. Overall, the dynamical evolution of the plasmoid instability appears to be similar to 
Cases D and E, including the temperatures reached at the main A-point (Eig. 7(b)), so that 
a saturation of the reconnection rate at similar levels is conceivable. 


The current sheet widths and aspect ratios for onset of the plasmoid instability are, 
approximately, 1200 m and 500 in Case D, 950 m and 730 in Case E, and 80 m and 1250 in 
Case E (see Eig. |^. The onset occurs at a higher aspect ratio when the sheet thins more 
dynamically. 


To summarize the runs for Model II without guide field, we find a basic analogy to 
the strong guide field cases in that the plasmoid instability is the dominant mechanism to 
realize fast reconnection, with rates similar to the strong guide field cases. Although our 
simulations confirm the strong thinning of the current sheet allowed by ambipolar diffusion 


and anticipated to yield an accelerated Sweet-Parker-like reconnection regime (Heitsch & 


Zweibel|2003a), this regime does not occur because the aspect ratio for onset of the plasmoid 


instability is reached by the end of the ambipolar diffusion-driven current sheet thinning. 







An accelerated Sweet-Parker-like reconnection due to radiative cooling of the current sheet 
does not occur in our model because we include a background heating term representing the 
chromospheric heating. 


3.3. Effects of the density stratification 

We hnd that the inclusion of gravity and density stratihcation in Model III has surpris¬ 
ingly little effect on the main results found for Models I and II, although the density varies by 
two orders of magnitude in the relevant height range y > Lq. Specifically, (1) the plasmoid 
cascading process shown in Fig. 2 is very similar to Models I and II; (2) the reconnection 
rate, current sheet width and temperature at the main X-point for Cases GM in Model III 
behave very similar to those for Cases A-C in Model I, respectively (Figures 5-7); (3) many 
fine slow-mode shock structures attached around the edges of the plasmoids appear, where 
the highest temperatures are produced; (4) the onset aspect ratios for Models I and III are 
the same, as shown in Fig. 8 and Table 1; (5) radiative cooling and ambipolar diffusion have 
little effect on the plasmoid instability in Model III; (6) the maximum out-flow velocity in 
Model III is the same as in Model I, about 40 kms“^. The similarity of the reconnection 
rate, which is normalized by the same value of as in Model I, and other quantities at 

the main X-point can be understood from the rather similar values of the Alfven velocity at 
that point in Models I and III, with the Alfven velocity initially being higher in Model III 
by only a factor ^1.2; this implies a similar dynamic behavior in the vicinity of the main 
X-point. 

There are some details in the reconnection process which are different in Models I 
and III. As shown in Figs. 8 and 9, the current sheet gradually becomes inclined to the 
side, especially after secondary islands appear. This is due to gravity acting on the dense 
plasmoids, which form slightly asymmetrically about the current sheet as a consequence of 
roundoff errors. Secondary islands in Model III start to appear later than in Model I for 
the same initial perturbation. In the period before secondary islands appear, the hottest 
plasma exists in the up-flow region in Model III (Figure 10), but in the down-flow region 
in Model I (not shown). Since low-density plasma can be heated more easily, the rapid 
decrease of the plasma density with increasing height in Model III is the main reason for this 
difference. For Case H, the hottest plasma usually appears around the edge of the plasmoids, 
but their center is still relatively much cooler (Figure 10), the reason being that the plasma 
density inside the plasmoids is higher and, therefore, radiative cooling is stronger. During 
the reconnection process, the main X-point in Model III gradually moves to a higher place 
compared to Model I, as seen in Figures 9 and 12. With radiative cooling, this results in a 
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higher temperature for Model III 80000 K in Case H vs. ~ 70000 K in Case B). 


3.4. Comparison with observations 


Micro-jets have been frequently observed in the low atmosphere (e.g., Liu et al. 2009; 
Morton 2012; Singh et al. 2012; Bharti, Hirzberger, & Solanki 2013). Magnetic reconnection 
is being considered as one of the mechanisms that may cause the ejection of these jets. Their 
velocity is usually found in the range of 10-150 kms“h The average speed of chromospheric 
jets around the edges of sunspots is around 30 kms“^ (Morton 2012). Fig. 12 shows that 
the upward outflow velocity Vy can reach about 40 kms“^ in the strong guide field model 
with radiative cooling and about 25 kms“^ in the zero guide field model. These upward 
reconnection outflow velocities are close to the typical velocities of chromospheric jets. Higher 
outflow velocities (by a factor 

/^O ^ 3) from a chromospheric current sheet with zero 

or weak guide held can be expected if the sheet is not overdense as in our Model II. The 
outflows become intermittent in the course of the plasmoid instability as relatively large 
merged islands are convected out of the sheet (Fig. [I^. This corresponds to observations of 
plasmoids in chromospheric jets (e.g., Singh et al.||2012 ). 

Many observations demonstrate that chromospheric jets can be heated to the transition 


region temperature (e.g., Teriaca, Curdt, & Solanki 2008). Morton (2012) pointed out that 
“slower” (v ~ 30 kms“^) jet events near sunspots involve 10^ K plasma and could play a role 
in heating the chromosphere and corona. Although radiative cooling is included in Cases B 
and E, these simulations yield temperatures in the current sheet region of several 10^ K 


(Figs. 3, 7, and 10), and the plasma immediately downstream of the slow-mode shocks and 
in turbulent regions between secondary islands is heated to temperatures of the lower transi¬ 
tion region, ~ 8 x 10^ K. Therefore, our simulations demonstrate that magnetic reconnection 
driven by the plasmoid instability in the low atmosphere contributes to chromospheric heat¬ 
ing. 


4. Conclusions and discussion 

Magnetic reconnection in the middle chromosphere, including radiative cooling, heating, 
and ambipolar diffusion effects, has been studied in the 2.5-dimensional MHD approximation. 
A Harris current sheet model with either a strong guide held (Models I and HI) or vanishing 
guide held (Model H) is considered. We have assumed that the plasma consisting of ions, 
electrons and neutral hydrogen atoms is strongly coupled and in ionization equilibrium state, 
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that the pressure is isotropic, and that the heat flux can be neglected. The ionization and 
radiative loss model by Gan fc Fang| ( |1990 ) is used, resulting in an initial ionization degree in 
the current sheet in the range ~ (0.2-0.5) %. A temperature-dependent magnetic diffusion 
coefficient adapted to the values in Khomenko & Collados (2012) is employed, leading to 
Lundquist numbers of S' ~ 10®-10^. 

The main numerical results can be summarized as follows: 


1. The plasmoid instability mediates a fast reconnection regime under chromospheric 
conditions for vanishing as well as for strong guide held. Reconnection rates of order 0.01- 
0.03, a factor ~ 20 above the Sweet-Parker rate and within the range of observationally 
inferred rates, are reached in either case, based on the classical Spitzer resistivity. The aspect 
ratio at the onset of the plasmoid instability at the supercritical values of the Lundquist 
number considered here is found to lie in the range 500-1250. 

2. AMR allowed us to resolve three levels of secondary islands and current sheet widths 
down to 5x ~ 30 m, where the one-fluid approximation is still valid. 


3. Ambipolar diffusion and radiative cooling have a significant influence on the re¬ 
connection process only in the model with vanishing guide field. Ambipolar diffusion then 
leads to the expected strong current sheet thinning (Brandenburg & Zweibel 1994), but 
the current sheet thins very fast, reaching the onset aspect ratio, such that an accelerated 
Sweet-Par ker reconnection regime does not develop, but rather the plasmoid instability sets 
in. The expected strong current sheet thinning in a Sweet-Parker regime due to radiative 
cooling ( Uzdensky fc McKinney||2011 ) does not occur in our simulations because we include 
a background heating term. This thinning can generally not be expected to be strong in 
the middle and lower chromosphere, where the temperature is hardly higher than twice the 
solar temperature minimum. Both ambipolar diffusion and radiative losses have an effect 
on the plasmoid instability for vanishing guide field, implying, respectively, a much earlier 
onset and a faster development, but apparently no significant increase in the reconnection 
rate. 


4. Many slow mode shocks are generated between secondary plasmoids and secondary 
fragments of the current sheet during the unstable reconnection process. Downstream the 
shock front regions (in the plasmoids), the temperature is highest and reaches T ~ 8 x 10^ K, 
i.e., lower transition region temperatures. This supports conjectures in the literature that 
fast magnetic reconnection is a candidate mechanism for the chromospheric heating. 

5. Upward outflow velocities in the course of the plasmoid instability reach 40 kms“^ 
and ~ 25 kms“^ in the strong and zero guide field cases, respectively. Outflow velocities 
higher by a factor ~ 3 are expected in the zero guide field case if the initial current sheet were 
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not chosen to be over dense. These values lie in the range observed for chromospheric jets. 
Dynamic plasmoids yield variations in the outflow, consistent with observations of plasmoids 
in chromospheric jets. 


Recent work by Leake et al. ( 2012 , 2013) has studied the effects of recombination on 
reconnection in the middle chromosphere in a two-fluid model (weak coupling between ions 
and neutrals). This also included the occurrence of the plasmoid instability. Reconnection 
rates of order 0.1 were found to result from recombination effects, while the plasmoid insta¬ 
bility gave only a minor additional increase of ~ 15%. Recombination within the plasmoids 
reduces their current, which may slow down the dynamics of the plasmoid instability. It 
thus appears that recombination must ultimately be included for a complete description of 
reconnection in partially ionized plasma. However, the small increase of the reconnection 
rate due to the plasmoid instability in Leake et al. ( 2012 ) results in part from their choice 
of a high magnetic diffusivity, 77 = 6 x 10 ^- 2 .4 x 10 ® m^s“^, which is about two orders of 
magnitude higher than a typical magnetic diffusivity in the middle chromosphere and did 
not change with temperature. Our diffusivity is r^o = 1.8 x 10^ m^ s“^ initially and decreases 
to ~ 2 X 10^ m^ s“^ as the current sheet is heated in the course of the plasmoid instability. 
This gives a difference in 77 of up to three orders of magnitude. In Leake et al. (2013) the 
diffusivity was computed for chromospheric temperatures and of order 3 x 10^ m^ s“^ before 
the onset of the plasmoid instability, much closer to our values. The reconnection rate then 
reached ~ 0.05, only a factor ~ 2 above our value. Thus, the high reconnection rate found 
in Leake et al. (2012) results in part from the high diffusivity chosen in this paper. The 


reconnection rate of 0.05 in Leake et al. (2013) is due to recombination effects and higher 
but comparable to the reconnection rate we And during the plasmoid instability. After the 
onset of the plasmoid instability in Leake et al. (2013), the reconnection rate shows a steep 
increase, which could only be followed up to ~ 0.06, since the current sheet width then 
reached the resolution limit. The steep increase can be taken as an indication that the 
plasmoid instability is potentially very important in the weak coupling description of chro¬ 
mospheric reconnection as well. Additionally we note that Leake et al.| used the magnetic 
held and Alfven velocity in the upstream region not very far from the current sheet to nor¬ 
malize the reconnection rate, while we used the asymptotic inflow values. Our reconnection 
rates increase by a factor 2-3 if we use the magnetic held and Alfven velocity at distances 
in the range 5-20 current sheet half widths ((2.5-10)5x). Since the initial plasma density at 
the initial main reconnection X-point {y ~ 1.5Lo) in our models (~ in Models I and 

III) is more than 10 times higher than that in their models (6 x 10^®m“^), and the initial 
temperature (To = 7000 K) is lower than in their models (Tq = 8500 K), our models represent 
the conditions at somewhat smaller heights in the chromosphere. 


Our simulations do not include thermal conduction, which can generally be expected 
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to be energetically less important than radiative cooling at chromospheric temperatures of 
~ 10^ K which dominate in the major part of our simulated volumes (Lin et ah 1992). 


Much higher temperatures and steep temperature gradients are built up in the plasmoids 
and at the slow-mode shocks and X-points. Conductive energy transport from these regions 
into the inflow regions and into the big plasmoid at the bottom end of the current sheet 
is largely directed across the magnetic held (conduction in the third direction, parallel to 
the held, is not relevant in our 2.5-dimensional simulations). Thermal conduction by the 
charged particles is very strongly suppressed in the direction perpendicular to the held; 
however, the contribution by the neutral component is not inhuenced. Deviations from local 
thermodynamic equilibrium (LTE) allow a neutral component to be present even at the high 
temperatures that we hnd in the plasmoids. To estimate the potential magnitude of thermal 


conduction for the conditions of our simulations, we follow Orrall & Zirker (1961). They hnd 


that the thermal conductivity across the held is provided by the neutral component in the 
whole temperature range relevant here and give the coefficient of thermal conductivity in a 
pure hydrogen gas as 

kBT 
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where Mh = 1 is the molecular weight and = 3.21 x 10“^° m^ is the neutral-neutral 
collisional cross section (note that Orrall & Zirker use to designate the cross section). 


Since the ionization degree is small and does not change with time, n^/{n^ -|- Ui) ~ 1 in our 
models. Under these conditions, the coefficient of thermal conductivity is evaluated as 


K 


2.59 X 


where T is measured in Kelvin and the dimension of/visJK ^m ^s 


(30) 


The thermal con¬ 


duction is given by 


Fc = V ■ (kVT). 


( 31 ) 


Our temperature gradient scale (slow-mode shock width) is typically ~ 500 m. For the tem¬ 
perature as high as T = 80000 K, the resulting thermal conduction Fc ~ 2.59 X 10-^T^/^S-^ 


is about 2.3 J m 


-3 C.-1 


This falls in the range of radiative cooling relevant in our models. 


Considering Model III, which encompasses densities uh ~ lO^^-lO^^ m“^ and corresponding 
ionization degrees Y\ ~ 10“^-10“^, and using a ~ 0.01 as a conservative estimate, plasma at 
80 000 K cools radiatively at rates of 0.35-350 J m“^ s“^. This suggests that neutral thermal 
conduction could limit the rise of the temperature in Model I and in the upper part of the 
volume in Model III, where the densities are relatively low. However, if the evolution of the 
ionization degree were also included, neutral thermal conduction can be expected to be neg¬ 


ligible. At T = 80000 K the ionization degree is very high, e.g., Orrall & Zirker (1961) give 
Un/ni ~ 3 X 10“® in their Table 5, so that the estimate of the neutral conductive heat flux 
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drops by this factor and is thus lower that the radiative losses by at least four orders of mag¬ 
nitude throughout the range of parameters studied in this paper. Additionally, more recent 
considerations of the collisional cross sections ( Khomenko fc Collados||2012 ) indicate a value 
for (Tnn higher by an order of magnitude than given by |Qrrall fc ZirkCT , The above estimates 
should be checked in future simulations of chromospheric reconnection when non-equilibrium 
ionization is also taken into account. 

Effects of potential relevance but not yet included in our model are radiative transport 
and non-LTE effects, non-equilibrium ionization effects, and at very small scales the Hall 
effect. Additionally, the complex topology of the chromospheric magnetic field has not yet 
been addressed. It is currently impossible to include all of them in a single model. Perhaps 
the most relevant steps beyond our approximations are the inclusion of recombination effects 
in a two-fluid description, more complex magnetic topologies guided by the observations, and 
radiative transport. Regarding the latter extension, we note that several estimates of the 


reconnection rate in plasmas with radiative cooling in Uzdensky & McKinney (2011) were 


found to not depend sensitively on the specific properties of the cooling process. Also, the 
effect of radiative losses is only a weak to moderate one in the computations presented here. 
These facts may be taken as an indication that the inclusion of radiative transport in the 
optically thick case would not strongly influence most of our results. 
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Table 1: Important parameters of the current sheet region in Models I, II and III. Initial 
values of Tq - temperature, Sq - Lundquist number, uho “ hydrogen density, riio - ion density, 
uao “ Alfven velocity, dio - ion inertial length, Anio ^ neutral-ion collision mean free path, 
and current sheet width (ions measured just before secondary islands appear. The height 
dependence of these parameters is significant only in Model III. 



ro(K) 

Ao 

nHo(m“3) 

nio(m“^) 

UAo(ni/s) 

dio(m) 

Anio(ni) 

^ons (^) 

Model I 

7000 

1.88 X 10® 

4.1 X 10^9 

2.3 X 10^^ 

3.32 X 10^ 

0.46 

4.68 

400 

Model II 

7000 

5.66 X 10® 

4.5 X lO^o 

7.5 X 10^^ 

1.00 X 10^ 

0.24 

0.88 

80-1200 

Model III 

{y = Lo) 

7000 

3.64 X 10® 

9.04 X lO^o 

1.05 X 10^® 

7.28 X 10^ 

0.22 

1.19 

400 

{y = l.SLo) 

7000 

1.18 X 10® 

8.59 X 10^*^ 

3.23 X 10^^ 

2.36 X 10^ 

0.39 

3.88 

400 

{y = 2Lo) 

7000 

3.83 X 10® 

8.17 X 10^8 

9.06 X 10^® 

7.66 X 10^ 

0.72 

12.58 

400 



(a)Model I, case B, t=9.9s 
Jz(A/nn2) 



Fig. 1.— (a) Field lines and current density (background color) in the whole simulation 
domain for the strong guide field model I with radiative cooling (Case B) at t = 9.9 s. (b) 
Field lines and plasma density p (background color) in the domain we focus on for Model III 
with gravity at t = 0 (with the initial perturbation included). 
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t = 41 .56s 


Fig. 2.— Field lines and current density showing four levels of plasmoids and fragmentary 
current sheet sections during the cascade to small scales for Model III, Case I at t = 41.56 s. 
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Jz(A/m^) T(K) 
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36231 
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1.694 


1.692 
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0.497 0.502 0.497 0.502 

x(l-o) x(l-o) 


Fig. 3.— Left: Field lines and cnrrent density of a plasmoid in Model I, Case B at t = 31.3 s. 
Right: Heating at the slow-mode shocks at the edge of the plasmoid. 
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x(Lo) 


Fig. 4.— Magnetic field components parallel and perpendicnlar to the slow-mode shock and 
current density along a cnt line at y = 1.691 Lq in Fig. (with a scale factor a — 0.00025 
adjnsting it to the range of the plot). 



t(s) 


l(s) 


t(s) 


Fig. 5.— Normalized reconnection rate in (a) Model I with strong guide field and without 
gravity, (b) Model II without guide field and without gravity, and (c) Model III with strong 
guide field and gravity. 
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(a) 



(b) 



t(s) 



t(s) 


Fig. 6.— Current sheet width (FWHM) at the main X-point in (a) Model I, (b) Model II 
and (c) Model III. 





Fig. 7 


Temperature at the main X-point in (a) Model I, (b) Model II and (c) Model III 
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Model I, cose B Model II, cose D Model II, cose E Model II, cose F Model III, cose I 
J,(A/m^) J,(A/m^) J,(A/nn^) J,(A/nn^) J,(A/m^) 


- 45.7 0.0 45.7 - 13.2 0.0 13.2 - 16.1 0.0 16.1 - 28.1 0.0 28.1 - 46.7 0.0 46.7 



xCLo) 


x(Lo) 


x(Lo) 


x(Lo) 


x(Lo) 


Fig. 8.— Field lines and current density at the times of plasmoid instability onset in Model I 
for Case B (t = 26.4 s), Model II for Cases D-F (t = 92.5 s, 84.0 s and 50.5 s, respectively), 
and Model III for Case I (t = 33.4 s), indicating the aspect ratio at the onset of the plasmoid 
instability in each case. 
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(a)Model II, case D 
t=53.5s J^(A/m^) t= 101.8s J^(A/m^) 


(b)Model II, case E 
t=53.5s J^(A/m^) t= 101.8s J^(A/m^) 


(c)Model II, case F 
t=53.5s J^(A/m^) t=101.8s J^(A/m^) 



(d)l\/1odel III, case G 
t=24.9s J 2 (A/m^) t=46.1s 


(e)Model III, case H 
t=24.9s t=46.1s J 2 (A/m^) 


(f)Model III, case I 
t=24.9s J 2 (A/m^) t=46.1s J^(A/m^) 


- 10.9 0 . 10.9 - 186.5 0 . 186.5 - 11.3 0 . 11.3 - 203.2 0 . 203.2 - 10.8 0 . 10.8 - 177.4 0 . 177.4 





x(Lo) 


x(Lo) 


x(Lo) 


x(Lo) 


x(Lo) 


x(Lo) 


Fig. 9.— Field lines and current density at two times in Model II for Cases D-F, and 
Model III for Cases G-I. The assignment of color to current density is adapted in each 
panel. The display is expanded in x direction to show the plasmoids clearly. 
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(a)Model II, case D 

t=55.5s T(K) t=101.8s T(K) 


(b)Model II, case E 
t=55.5s T(K) t=101.8s T(K) 


(c)Model II, case F 
t=55.5s T(K) t=101.8s T(K) 


4464 6932 9400 3604 6495.5 9387 

2.0 E 





0.485 0.500 0.515 0.485 0.500 0.515 0.485 0.500 0.515 0.485 0.500 0.515 0.485 0.500 0.515 0.485 0.500 0.515 


x(Lo) x(Lo) 

(d)Model III, case G 

t=24.9s T(K) t=46.1s T(K) 


>^(Lo) x(Lo) 

(e)Model III, case H 

t=24.9s T(K) t=46.1s T(K) 


x(Lo) x(Lo) 

(f)Model III, case I 

t=24.9s T(K) t=46.1s T(K) 


6810 13106 19402 4621 47026 89431 6514 9957.5 13401 4518 45573 86628 6807 13103.5 19400 5176 48605.5 92035 




0.485 0.500 0.515 0.485 0.500 0.515 

x(Lo) x(Lo) 


0.485 0.500 0.515 0.485 0.500 0.515 0.485 0.500 0.515 0.485 0.500 0.515 

x(Lo) x(Lo) x(Lo) x(Lo) 


Fig. 10.— Field lines and temperatnre for the same cases and times as in Fig. The 
assignment of color to temperature is also adapted in each panel. 
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(a) (b) 

t=16.6s t = 27.1s 




t=27.1s, E*oz=[- 28.45, 1.69]mT/s 


Fig. 11.— Profiles in x direction of the 2 ; component of the ambipolar diffusion field Eadz 
(solid black line) and — (jyV x B)^ (red dash line) for Model I, Case C at the main X-point at 
(a) t = 16.6 s and (b) t = 27.1 s. (c) Spatial distribution of Eadz around the main X-point 
Sit t — 27.1 s. 
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Fig. 12.— Vertical velocity component for Model III, Case H. 








